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ABSTRACT 

New and upgraded radio interferometers produce data at massive rates and will require significant 
improvements in analysis techniques to reach their promised levels of performance in a routine manner. 
Until these techniques are fully developed, productivity and accessibility in scientific programming 
environments will be key bottlenecks in the pipeline leading from data-taking to research results. 
We present an open-source software package, miriad- python, that allows access to the MIRIAD 
interferometric reduction system in the Python programming language. The modular design of 
MIRIAD and the high productivity and accessibility of Python provide an excellent foundation for 
rapid development of interferometric software. Several other projects with similar goals exist and we 
describe them and compare miriad- python to them in detail. Along with an overview of the package 
design, we present sample code and applications, including the detection of millisecond astrophysical 
transients, determination and application of nonstandard calibration parameters, interactive data 
visualization, and a reduction pipeline using a directed acyclic graph dependency model analogous 
to that of the traditional Unix tool nnake. The key aspects of the nniriad-python software project are 
documented. We find that nniriad-python provides an extremely effective environment for prototyping 
new interferometric software, though certain existing packages provide far more infrastructure for some 
applications. While equivalent software written in compiled languages can be much faster than Python, 
there are many situations in which execution time is profitably exchanged for speed of development, 
code readability, accessibility to nonexpert programmers, quick interlinking with foreign software 
packages, and other virtues of the Python language. 
Subject headings: Data Analysis and Techniques 



1. INTRODUCTION 

Advances in the fields of digital computing and com- 
mercial wireless communication have fueled an explosion 
of innovation in radio interferometry. New and upgraded 
facilities such as the Allen Telescope Array (ATA; Welch 
et al. 2009), the Low- Frequency Array (Kassim et al. 
2004), the Precision Array for Probing the Epoch of 
Reionization (PAPER; Parsons et al. 2010), the Murchi- 
son Wide-Field Array (MWA; Lonsdale et al. 2009), the 
Karl G. Jansky Very Large Array (formerly EVLA; Per ley 
et al. 2011), the Westerbork Synthesis Radio Telescope 
(WSRT) Apertif project (Verheijen et al. 2008), the Aus- 
tralian Square Kilometer Array Pathfinder (DeBoer et al. 
2009), and MeerKAT (Jonas 2009) have sophisticated 
designs including large-number-of-small-dishes architec- 
tures, multipixel or phased-array feeds, wide bandwidths, 
and phased array substations. Several of these facilities 
are pathfinders for the proposed Square Kilometer Array 
(SKA; CarilH & RawHngs 2004). They all aim to push the 
limits of interferometric techniques to image large fields 
of view with unprecedented spatial, spectral, temporal, 
and polarimetric fidelity, at unprecedented data rates. 

Achieving these goals will require substantial amounts 
of new software to process the data coming out of these fa- 
cilities. The data rates present a challenge in and of them- 
selves: SKA-class facilities will require exaflop-scale com- 
puting (Cornwell & Humphreys 2010). Another challenge 
is the development of the necessary new algorithms, which 
is certain to require extensive experimentation. Tech- 
niques already under investigation include K;-projection 
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(Cornwell et al. 2008), A-projection (Bhatnagar et al. 
2008), multi-scale mult i- frequency CLEAN (Rau & Corn- 
well 2011), delay/delay-rate filtering (Parsons & Backer 
2009), space-alternating generalized expectation (SAGE) 
maximizing calibration (Kazemi et al. 2011), generalized 
measurement-equation-based instrumental modeling (No- 
ordam & Smirnov 2010), scale- invariant rank detection 
of radiofrequency interference (RFI; Offringa et al. 2012), 
subspace-tracking RFI mitigation (Ellingson & Hampson 
2002), visibility stacking (Hancock et al. 2011), bispec- 
tral pulse detection (Law & Bower 2011) and improved 
sourcefinding tools (Whiting 2012), to name a few. Until 
the next generation of algorithms is thoroughly explored, 
the quality of many results will be set not by the capabil- 
ities of observatory hardware but by the sophistication of 
the reduction pipeline that can be brought to bear before 
publication: software-limited science. 

Although improvements in software development effi- 
ciency are always desirable, they're particular salient now 
as the next generation of radio interferometers comes on- 
line. We discuss the efficiency of a programming environ- 
ment in terms of productivity^ which we take to measure 
the amount of useful functionality that a programmer 
can implement per unit time, and accessibility^ which 
measures how much effort it takes for non-experts to be- 
gin successfully working within the environment without 
hand-holding. Although many studies attempt to quan- 
tify these metrics (e.g., Petersen 2011, and references 
therein), our discussion will remain qualitative. These 
attributes are relevant in situations in which developer 
resources are constrained: a more productive environment 
allows individual developers to accomplish more, while a 
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more accessible environment has a lower barrier to entry 
from informal contributors. 

Interpreted, dynamic programming languages can pro- 
vide such an environment. While several of these exist, 
the language Python^ in particular has seen a broad up- 
take in the astronomical community over the past decade 
(e.g., Barrett & Bridgman 1999; Greenfield & White 2000; 
Blakeslee et al. 2003; Kettenis et al. 2006; Magee et al. 
2007). It is intended to be easy to learn and offers conve- 
niences including object-oriented programming, lambda 
expressions, exception handling, and an enormous soft- 
ware ecosystem, including interactive interpreters (e.g., 
IPython; Perez & Granger 2007), visualization tools (e.g., 
matplotlib; Hunter 2007), file format interfaces (e.g., py- 
FITS; Barrett & Bridgman 1999), database interfaces (e.g., 
SQLite^), statistics routines, web service support, and so 
on. The existence of this ecosystem and the rapid rise in 
the popularity of Python in the astronomical community 
in particular speak to its productivity and accessibility. 

The programming environments of "classic" astronomi- 
cal software packages tend to be intimidating and frus- 
trating in comparison. Two of the major traditional 
radio interferometric reduction packages, MIRIAD (Mul- 
tichannel Image Reconstruction, Analysis, and Display; 
Sault et al. 1995) and AIPS (Astronomical Image Pro- 
cessing System; Greisen 2002), are largely implemented in 
FORTRAN-77, and in both cases the process to go from 
source code to usable installed executable code is compli- 
cated and fragile. In light of the clear need for improved 
interferometric algorithms, it should be no surprise that 
a variety of projects have sought to build on or replace 
these systems with more modern ones. Perhaps the most 
prominent such undertaking is the Common Astronomy 
Software Apphcations (CAS A; McMulhn et al. 2007), the 
successor to AIPS. Other packages include Obit (Cotton 
2008b), AIPY (Astronomical Interferometry in Python^), 
MeqTrees (Noordam & Smirnov 2010), and several more 
narrowly-targeted binding layers (e.g., pyramid; Mehringer 
& Plante 2004). 

To this list, we add miriad- python, a creatively- named 
package exposing MIRIAD tasks and subroutines in 
Python. Its design and philosophy are discussed (§2) 
and compared those of related projects (§3), including 
the ones mentioned above. We then describe the imple- 
mentation of miriad-python (§4) and provide some very 
brief examples of its use (§5). Some applications in which 
it has been used are presented, including the detection of 
millisecond astrophysical transients, determination and 
application of nonstandard calibration parameters, inter- 
active data visualization, and pipeline processing (§6). 
We document the nature of miriad-python as a software 
project (§7), discuss some of its performance characteris- 
tics (§8), and finally summarize (§9). 

2. DESIGN CONSIDERATIONS IN miriad-python 

The miriad-python project is intended to allow conve- 
nient access to MIRIAD tasks, datasets, and subroutines. 
Although it aims to ease the rapid development of new in- 
terferometric algorithms, it does not provide any nontriv- 
ial algorithms itself. As such, the focus of miriad-python 
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is very narrow: it provides the best possible interfaces 
for access to MIRIAD infrastructure and is as fiexible as 
possible regarding what is done with it. To use terminol- 
ogy dating back to at least the design of the X Windows 
system, it provides mechanism but not policy (Scheifler 
& Gettys 1986). We believe that the wide range of ap- 
plications described in §6 is evidence of the power of this 
approach. 

Although work on miriad-python was initially inspired by 
a practical desire for more efficient scripting of MIRIAD 
reductions, the MIRIAD package is a good fit to the 
larger goals of the miriad-python project. MIRIAD offers 
its wide variety of tools and algorithms in a modular task 
architecture, and its data formats are simple and efficient. 
In particular, it is fairly straightforward to (ab)use the 
streaming MIRIAD u-v data format for novel applications 
(see, e.g., AIPY, §3.2). Other packages are built upon 
MIRIAD for similar reasons (e.g., Teuben 2011; Pound & 
Teuben 2012). MIRIAD is routinely used to process data 
from facilities including the ATA, the Berkeley-Illinois- 
Maryland Association array (BIMA; Welch et al. 1996), 
the Combined Array for Research in Millimeter/submil- 
limeter Astronomy (CARMA; Woody et al. 2004), the 
Australia Telescope Compact Array ( ATC A) , and the Sub- 
millimeter Array (SMA; Ho et al. 2004). MIRIAD is not a 
monolithic project: there are at least two nontrivially di- 
vergent codebases maintained for use with the ATCA and 
CARMA, with some sharing of modifications between 
the them, miriad-python is referenced to the CARMA 
MIRIAD codebase. Recent work on MIRIAD includes 
complete support for 64-bit file offsets and pointers, allow- 
ing imaging of arbitrarily large datasets; integration of the 
wcslib library (Calabretta 2011) for more comprehensive 
support of coordinate manipulations; and of course bug 
fixes, new features, and documentation improvements. 

When creating a package such as miriad-python that 
interfaces with a lower-level one, one must decide how 
closely to hew to the APIs (application programming 
interfaces) of the original package. In miriad-python, the 
APIs have been significantly reworked: they are object- 
oriented and aim to take full advantage of builtin language 
features. When they meet the goal of providing mech- 
anism rather than policy, substantial new features are 
implemented in the Python layer. One example is ap- 
proximate cryptographic hashing of u-v datasets, which 
can be used to quickly and fairly robustly check for mod- 
ifications to the data regardless of file modification times. 
(This is useful in dependency-tracking pipelines, cf. §6.4.) 
To write useful programs, miriad-python authors must be 
familiar with the semantics of MIRIAD datasets or tasks, 
and good knowledge of Python is helpful, but they do not 
need to understand the design of the MIRIAD subroutine 
library. 

Another consideration for authors is whether to link 
to a separate installation of the lower-level software or 
to compile and install their own version of it, bundling a 
copy of its source code with their own. We have chosen to 
go the former route with miriad-python. Although miriad- 
python offers very different APIs and new features, we 
concieve of it as fundamentally a layer above MIRIAD 
and not a standalone package: if the user has customized 
their MIRIAD installation in some way, for instance, it's 
appropriate for miriad-python to reffect that customization, 
and not override it with its own copy of the subroutine 
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hbrary. Another reason is that MIRIAD is stih evolving, 
and tracking changes from an upstream codebase to a 
forked copy can be tedious and error-prone. Another is 
that miriad-python is a general-purpose interface to not 
just the MIRIAD subroutine library but its task collection, 
so the entire, substantial, MIRIAD codebase would need 
to be bundled. Finally, we expect that miriad-python 
users will be existing MIRIAD users who are likely to 
already have installed MIRIAD on their systems, so the 
convenience of bundled source code is lessened. 

To touch on a broader philosophical issue, it may be 
argued that accessibility is not a desirable feature when 
it comes to interferometric software. Interferometric anal- 
ysis tends to be subtle, and most people write buggy 
code (see §7.5). On the other hand, scientific research is 
a fundamentally creative and exploratory process often 
requiring new or improved techniques, and that is clearly 
the case in the domain of interferometry at present. The 
tension here is a variation on the well-explored "cathe- 
dral versus bazaar" theme first described by Raymond 
(1999). Considering the amount of algorithmic devel- 
opment necessary to fully exploit next-generation radio 
interferometers, we take the stance that the empowerment 
of interferometric software users is a good thing. 

3. COMPARISONS TO RELATED PACKAGES 

As alluded to in §1, there are several other efforts un- 
derway to provide developer-friendly environments for 
producing interferometric software. In this section we 
attempt to situate miriad-python among them. We con- 
sider "bindings", which do not provide substantial new 
algorithms, and "high-level" packages, which do. 

3.1. Bindings 

Several other packages expose MIRIAD functionality 
in the Python language, but we found them unsatisfac- 
tory. The pyramid package (Mehringer & Plante 2004), 
providing a module named Miriad, includes similar func- 
tionality to that of mirexec in a somewhat less structured 
manner. The modules mirlib (via the WSRT) and miri ad- 
wrap both bind the MIRIAD subroutine library. Unlike 
miriad-python, they provide a fairly direct mapping to the 
MIRIAD subroutines and do not extensively "Pythonify" 
the API. These packages are less mature than miriad- 
python and do not appear to be under active development. 
None of them includes substantial documentation, high- 
level convenience features, or an integrated system for 
invoking both MIRIAD tasks and subroutines. 

There also exist several other packages that appear to 
share the same motivations as miriad-python but involve 
different technologies. Other dynamic languages can be 
used. For instance, MIRIAD functionality can be accessed 
in Ruby"^ with the MIRIAD-Ruby^ or mirdl^ packages. The 
two packages, written by the same author, cover the link- 
versus-bundle tradeoff discussed above: MIRIAD- Ruby 
includes a portion of the MIRIAD source code, while 
mirdi links to a separate installation of the full MIRIAD 
libraries. Although it is beyond the scope of this paper to 
compare the merits of the languages, we note that Ruby 
has seen less uptake in the astronomical community than 

^ http://www.ruby-lang.org/ 

^ http: //purl . org/net/pkgwpub/miriad-ruby 

^ http : //purl . org/net/pkgwpub/mirdl 



Python but is starting to make inroads (e.g., Lammers 
et al. 2002; Fuentes et al. 2007; Gutierrez-Kraybill et al. 
2010). 

Other packages build on top of different interferometry 
frameworks. ParselTongue (Kettenis et al. 2006) provides 
a Python interface to AIPS, building on top of parts 
of the Obit package (Cotton 2008b), described in the 
next subsection. The pyrap project^ provides a Python 
interface to the core C++ support libraries of CASA. 

In many cases, the choice of the appropriate binding 
will be constrained by external limitations in either the 
high-level language or the underlying interferometry pack- 
age. Our reasons for preferring Python and MIRIAD, 
respectively, are sketched in §1 and §2. We emphasize the 
simplicity and efficiency of MIRIAD 's data formats as 
being important factors in facilitating algorithmic exper- 
imentation. Another distinction between miriad-python 
and alternative bindings is its rate of development. Even 
in otherwise stable and mature packages, documentation 
improvements and subtle bugfixes tend to require fairly 
frequent modifications. In the calendar year of 2011, there 
were 5 commits to the mirdi version control system (VCS), 
10 commits to the pyrap VCS, and 110 commits to that of 
miriad-python. The last known update to ParselTongue 
was in 2010 and the last known update to pyramid was 
in 2008. Finally, we believe that miriad-python acquits 
itself well when it comes to the subjective factors of API 
cleanliness, code readability, and code quality. 

3.2. High-Level Packages 

The projects described in this subsection are more 
ambitious than miriad-python and include significant func- 
tionality along with Python-based environments: they 
provide policy, as well as mechanism (cf. §2). While 
the packages described in the previous subsection are 
intended to provide suitable bases for development for a 
wide range of applications, the ones described below may 
not be appropriate for certain applications, depending 
on both technical and architectural factors. All of these 
packages are actively developed. 

CASA (McMulhn et al. 2007) is a substantial pack- 
age including high-level Python interfaces to lower-level 
routines implemented in C++ and FORTRAN. The first 
iteration of CASA, known as AIPS++ (McMullin et al. 
2004), included flexibility, ease of programming, and mod- 
ifiability in its mission statement (AIPS++ Steering Com- 
mittee 1993). We find that these are difficult to achieve 
in practice because most CASA functionality is imple- 
mented in fairly large C++ modules, with the built-in 
Python interfaces providing mostly coarse-grained ac- 
cess to data and algorithms; the combination of the core 
CASA libraries with pyrap is more friendly to Python 
development. One minor inconvenience is that CASA 
bundles its own Python interpreter, so that site-specific 
packages accessible to the system interpreter sometimes 
need to be rebuilt and reinstalled in order to be accessible 
to CASA's version. CASA's imager is quite flexible and 
can combine multiple advanced algorithms, for instance 
simultaneous multi-scale, multi-frequency deconvolution 
(Rau & Cornweh 2011). 

Obit is billed as "a development environment for astro- 
nomical algorithms" (Cotton 2008b). Although many of 

^ http : //purl . org/net/pkgwpub/pyrap 
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its fundamental features are general-purpose, Obit's main 
emphasis is radio astronomy, and in particular it includes 
several algorithms important to low-frequency radio in- 
terferometry. The package includes a set of base libraries 
implemented in C, a set of tasks built on top of them, and 
a Python layer for interacting with both of these. This 
layer is based on ParselTongue, although the two projects 
remain separate (W. D. Cotton, priv. comm.). Two as- 
pects of the Obit design are of particular note. Firstly, 
Obit uses AIPS and FITS file formats for data storage, so 
it interoperates seamlessly with certain existing packages 
and developer effort need not be spent reimplementing 
existing tools. Secondly, the Obit libraries are architected 
to allow multithreaded execution, and this capability has 
been explored during Obit's development. In one example 
problem, the time taken to generate a complex image was 
reduced by a factor of 6.5 when the number of CPU cores 
used for computations was increased from one to twelve 
(Cotton & Perley 2010). 

The AIPY system represents a sharper break from 
tradition than those previously mentioned. It is a stan- 
dalone package implementing many advanced synthesis 
and calibration techniques in Python to meet the needs 
of the low-frequency dipole array PAPER. Its design em- 
phasizes flexibility and is Python-centric: the primary 
interfaces are the Python ones, and code from preexist- 
ing packages such as HEALPix (Gorski et al. 2005) and 
XEphem (Downey 2011) is liberally borrowed to provide 
fast, well-tested implementations of many features. AIPY 
can use the MIRIAD format for u-v data storage, and 
hence includes an internal Python binding of the relevant 
portions of the MIRIAD subroutine library. A particu- 
larly interesting strength of AIPY is its infrastructure for 
creating maps of the whole sky at once. 

The MeqTrees project (Noordam & Smirnov 2010) also 
uses Python to provide a productive environment for ex- 
ploring new interferometric algorithms using the measure- 
ment equation formalism (Hamaker et al. 1996; Smirnov 
2011). Its input/output (I/O) and astronomical infras- 
tructure are built on the core CAS A libraries and pyrap, 
with certain numerical routines implemented in C++. 
MeqTrees provides a substantial infrastructure for devel- 
oping and applying sophisticated calibration methodolo- 
gies with numerous visualization features. As an example 
of its power, MeqTrees has been used to generate a noise- 
limited radio astronomical image of the source 3C 147 
with a dynamic range of 1.6 x 10^ (Noordam & Smirnov 
2010). 

4. IMPLEMENTATION 

miriad- python provides three Python modules: mir- 
task, mirexec, and miriad. The first of these uses C and 
FORTRAN-77 based extension modules to interface with 
the MIRIAD subroutine library, providing most of the 
low- level functionality that makes miriad-python useful. 
Numerical functionality is provided by NumPy^ and link- 
ing with FORTRAN-77 is accomplished using f2py (Peter- 
son 2009). The second module, mirexec, is Python-only 
and provides a uniform framework for executing MIRIAD 
tasks inside the Python language. The third, miriad, is 
Python-only and provides a simple abstraction for refer- 
encing and manipulating MIRIAD datasets, with hooks 

^ http://numpy.scipy.org/ 



into mirtask and mirexec to provide convenient access 
to certain common operations. The three modules are 
loosely-coupled such that they do not actually depend on 
each other, although most nontrivial tasks end up using 
functionality from all three modules. The relationships 
of the APIs are diagrammed in Figure 1. 

As mentioned above, the MIRIAD APIs are significantly 
reworked and expanded. In most cases. Python's runtime 
type information is used to automatically choose data 
types when performing lowlevel I/O. Because these types 
are important for task interoperability, they must be 
specified explicitly when writing new data. MIRIAD 's 
internal errors are mapped into Python exceptions using 
a mechanism based on the C library functions setjmp and 
longjmp. 

miriad-python also provides a program that can serve 
as a help viewer for either existing MIRIAD tasks or 
programs written in Python, in the latter case using 
standard Python "documentation strings" written in the 
standard MIRIAD documentation format. Tasks written 
in Python can thus seamlessly merge with traditional 
compiled MIRIAD tasks from the user standpoint. 

5. EXAMPLE CODE 

Although it is impractical to provide extensive code 
examples in this paper, in this section we give a few 
brief ones. These are necessarily simplistic and do not 
demonstrate the most interesting capabilities of miriad- 
python. As discussed below (§7.4), the miriad-python 
source distribution contains longer samples^. We assume 
a familiarity with typical MIRIAD usage and the Python 
language. Many valuable coding practices are ignored 
below for the sake of concision. 

MIRIAD tasks can be executed with the mirexec module. 
The following example demonstrates the parallel execu- 
tion of multiple tasks: the launch method starts a task 
process, returning a handle to it (a subclass of Python's 
subprocess.Popen), and the handle's wait method pauses 
the caller until the process completes. (For I/0-intensive 
operations as shown here, it may actually be faster to 
execute the tasks serially to avoid disk thrashing.) 

import sys 

from mirexec import TaskUVAver 

prochandles = [ ] 

for p in sys.argv[l:]: 

taskdesc = TaskUVAver (vis=p, out=p-|-".avg", 

interval=10) 
prochandle = taskdesc. launch () 
prochandles. append (prochandle) 

for prochandle in prochandles: 
prochandle. wait () 



^ Currently browseable online at https://github.com/pkgw/ 
miriad- python/tree/master/examples. 
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The above example is intentionally verbose. To stim- 
ulate the reader's imagination, we provide the following 
variation that uses Python's "list comprehension" syntax 
and avoids the use of several local variables: 

import sys 

from mirexec import TaskUVAver 

for ph in [TaskUVAver (vis=p, out=p-|-".avg", interval=10) 
.launch () for p in sys.argv[l:]]: 
ph. wait 

The following example demonstrates how the "header 
variables" of MIRIAD datasets may be accessed using 
miriad- python. The variable vispath is some string giv- 
ing the filesystem path of a MIRIAD u-v dataset. As 
mentioned in §4, storage types of variables are detected 
automatically on read, but must be specified explicitly 
on write to enforce consistency with preexisting MIRIAD 
tasks. 

import miriad, numpy as np 

vispath = ... 

ref = miriad. VisData (vispath) 

handle = ref.open ("rw") 

ncorr = handle. getScalarltem ("ncorr") 

handle. setScalarltem ("mine", np. double, np.pi * ncorr) 

The following example replaces the contents of an image 
with the magnitudes of its fast Fourier transform (FFT). 
This can be used, for instance, to recover the weights 
used in the imaging process from the dirty beam image. 
The variable impath is analogous to vispath. The value 
of the variable whichplane indicates that the first image 
plane should be selected. 

import miriad 

from numpy import abs, float32 

from numpy.fft import fFt2, fftshift, ifftshift 

impath = ... 

whichplane = [ ] 

handle = miriad. ImData (impath). open ("rw") 

d = handle. readPlane (whichplane). squeeze () 

d = abs (ifftshift (fFt2 (fftshift (d)))) 

handle. writePlane (d.astype (float32), whichplane) 

handle. close () 

Finally, this example inserts information about the 
synthesized beam shape of an image into a preexisting 
relational database. 

import miriad, sqliteS 

impath = ... 

handle = miriad. ImData (impath). open ("rw") 

dbpath = ... 

db = sqliteS. connect (dbpath) 

info = [impath] 

for item in "bmaj bmin bpa". split (): 

info. append (handle. getScalarltem (item)) 

db.execute ("INSERT INTO data VALUES (?,?,?,?)", info) 
db. commit () 
db. close 

Such functionality would be difficult to achieve using 
FORTRAN. We emphasize that in the interests of sim- 
plicity and clarity, we have avoided some of the most 



important features of the Python language such as object 
orientation, first-class functions, and exception handling. 

6. APPLICATIONS 

The fast Python development cycle makes it relatively 
easy to produce new algorithms and tools for astronomical 
data analysis. In this section we present some example 
applications of miriad- python. 

6.1. Algorithms for Millisecond Transients 

We have investigated with new algorithms to study 
millisecond transients in the visibility domain using data 
from the ATA and VLA (Law et al. 2011; Law & Bower 
2011). Since interferometers have not traditionally oper- 
ated at this time scale, many interferometric packages lack 
features needed to work in this regime, e.g., the ability 
to dedisperse visibilities and high-resolution timestamps. 
We used miriad- python to develop new ways to visualize 
and search these unusual data streams for transient radio 
sources. 

Figure 2 demonstrates imaging of millisecond pulses 
from the Crab pulsar (B0531+21) using a Python-based 
toolchain. We observed the Crab pulsar at a time resolu- 
tion of 1.2 ms and frequencies between 720 to 800 MHz 
using the ATA (Law et al. 2011). At this time resolu- 
tion, the dispersion introduced by the interstellar medium 
(ISM) delays the pulse arrival time by a few tens of mil- 
liseconds across the band. The upper panel of Figure 2 
shows this effect in a spectrogram formed by summing 
visibilities, effectively forming a synthesized beam toward 
the Crab pulsar. The spectrogram shows that the dis- 
persive delay follows a quadratic shape that has a slope 
consistent with that expected from the Crab pulsar. After 
identifying the time and dispersion of the pulse, we can 
image it to see if it is consistent with a point source in the 
location of the Crab pulsar. The u-v input/output facili- 
ties of miriad-python are used to write a new dedispersed 
visibility dataset. The lower panel of Figure 2 shows 
an image made from visibilities during the pulse. The 
dispersive delay is large enough that it must be corrected 
in order to detect the source. 

Low-level access to visibility data in miriad-python also 
makes it possible to experiment with new algorithms. 
We've used this capability to develop a technique to de- 
tect millisecond transients based on an interferometric 
closure quantity called the bispectrum (Cornwell 1987; 
Law & Bower 2011). The bispectrum is formed by multi- 
plying three visibilities from baselines that form a closed 
loop. This product is sensitive to transients anywhere 
in the field of view, which makes it powerful for surveys. 
However, most software packages only use the bispectrum 
for calibration, so none have low-level access to functions 
for prototyping algorithms. Figure 3 shows millisecond 
light curves toward pulsar B0329+54 made with miriad- 
python. With direct access to dedispersed visibilities, we 
were able to compare traditional beamforming with the 
bispectrum technique. In this case, we can detect pulses 
for four of the five rotations of the pulsar during this 
observation. Being able to compare the techniques on 
the same pulses allowed us to show that the bispectrum 
responds more strongly than beamforming, confirming an 
unusual theoretical property of the bispectrum (Law & 
Bower 2011). 
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6.2. Retroactive SEFD Calibration 

The ATA lacks online measurement of the system tem- 
perature (Tgys) or system equivalent flux density (SEFD), 
information that is important for assessing data quality 
and thermal noise limits. We have written a task, ca I noise, 
to assess per-antenna SEFDs after-the-fact using obser- 
vations of a bandpass calibrator and the assumption that 
the variance across the channels of each spectral window 
is entirely thermal. 

Each ATA antenna has two orthogonally linearly po- 
larized feeds. The signal from each of these is identified 
as originating from a given antenna-polarization pair- 
ing ("antpol"). The SEFDs of the two feeds on a single 
antenna may differ appreciably, and all of the SEFDs 
may be time- variable. Rather than directly computing 
SEFDs, ca I noise determines calibration coefficients relat- 
ing SEFDs to the RMS value across the spectral window 
of each antpol's uncalibrated autocorrelation amplitude, 
a value we refer to as the RARA (raw autocorrelation 
RMS amplitude). While the value of this coefficient is 
assumed to be static over the course of an observing ses- 
sion, the derived SEFD value will vary proportionally 
with the RARA. Although the RARA is a function of 
both system noise and a gain factor, the time variation 
in the latter is small enough so that the derived SEFDs 
will be sufficiently accurate. 

ca I noise takes as an input a dataset with a gains (and op- 
tionally bandpass) table. The dataset is scanned through 
once, without applying the gains table, to accumulate 
RARA values for each timestamp and antpol. The dataset 
is then reread, with the gains being applied, to compute 
per-baseline SEFDs: 

1 



SEFD = - (cr^ + cFi) r]V2Aur, 

where a^ and a^ are the standard deviation in the real and 
imaginary parts, respectively, of the calibrated visibilities 
across the spectral window (thus in Jansky units), r] 
is a tabulated efficiency of the correlation (depending, 
e.g., on the number of bits used in digitization), Av is 
the channel width in the spectral window, and r is the 
integration time of the record. Per-antenna calibration 
coefficients are then calculated using a least-squares fit 
of the baseline-based values assuming a simple geometric 
dependence: 



SEFD^,^,^ = VqMM^7^c^R^ARA~I, 

where the q are the desired calibration parameters. Out- 
lier values (possibly due to RFI or hardware failures) are 
identified by user-specified chpping limits and iteratively 
removed from the fits. MIRIAD's on-the-fly calibration 
system does not support SEFD information, so the q are 
then recorded in a textual table on disk; a companion task 
reads the table along with an input dataset and creates a 
new dataset with the SEFD information inserted. (For 
technical reasons, a constant Tgys is assumed and varying 
values of the MIRIAD variable jyperk are inserted into 
the dataset, where SEFD = Tgys -jyperk and thus jyperk 
is related to the effective area of each receiving element.) 
calnoise is written in Python and uses nniriad-python 
to read the visibility data and for miscellaneous astro- 
nomical routines. It uses NumPy for its numerics and 
other Python libraries for optional plotting of diagnostic 



information, least-squares solving, storage of numerical 
metadata, and integration into the ATA pipeline described 
below (§6.4). The main implementation comprises ~600 
statement lines of code (i.e., excluding whitespace and 
comments). The chief benefits of using Python and miriad- 
python were quick creation of the task skeleton (thanks 
mainly to high-level data structures), rapid turnaround 
during refinement (thanks mainly to the lack of a compi- 
lation step) , and easy investigation of algorithmic tweaks 
(thanks mainly to easy access to other libraries, e.g., only 
one line of code needed to interactively plot variables) . 

6.3. Interactive u-v Data Visualization 

Data visualization is an important part of the develop- 
ment of any processing pipeline. This is especially true for 
radio interferometers, in which the data undergo complex 
transformations, are high-dimensioned, and closed-loop 
instrumental modeling is often a key aspect of the reduc- 
tion. We have used miriad-python to implement a powerful 
u-v data visualizer, first described in Williams (2010). As 
shown in Fig 4, the main display is a dynamic spectrum 
of the visibilities on a baseline as a function of frequency 
and time. The user can quickly switch between different 
displays (real, imaginary, amplitude, phase), apply var- 
ious processing steps on-the-fly (e.g., average, rephase), 
and navigate through the dataset. There is also substan- 
tial support for visualizing and creating data flags when 
necessary, although it is widely recognized that manual 
flagging of data is rapidly becoming an impossible task 
as data rates increase (Keating et al. 2010). 

The use of Python is a key aspect to the u-v visual- 
izer because it allows easy combination of the existing 
MIRIAD libraries with very different software, in this 
case the modern graphical toolkits Cairo^^ and GTK+-'^-'^. 
Cairo provides routines for fast rendering of the grid- 
ded visibility data, while GTK+ provides a higher-level 
widget library that allows new user interactions to be 
implemented quickly. Achieving this functionality with- 
out the use of these (or similar) libraries would be a 
massive undertaking: Cairo and GTK+ comprise 240000 
and 560000 lines, respectively, of well-tested, efficient C 
code. The main visualizer codebase, on the other hand, 
comprises only 2900 statement lines of Python. 

6.4. An ATA Reduction Pipeline 

Modern and next-generation radio observatories are 
expected to produce data at rates significantly higher 
than those of older facilities (e.g., Cornweh & Humphreys 
2010). Under these conditions, rapid automated pro- 
cessing of data is changing from a luxury to a necessity 
(Keating et al. 2010). Python is often identified as a good 
system in which to construct pipelines to perform such 
processing, thanks to its amenability to implementing 
high-level application logic, revisability, and ability to 
glue together existing tools and libraries (Sanner 1999; 
Myers et al. 2007; Perez et al. 2011). We have used 
nniriad-python to implement a commensal observing sys- 
tem (Wilhams 2012) for the ATA Galactic Lightcurve 
and Transient Experiment (AGILITE; Williams et al., 
2012, in prep.) and just such a pipeline for processing 

^^ http://cairographics.org/ 
11 http://gtk.org/ 
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the data. In the latter case, the mirexec module is used 
extensively to invoke existing MIRIAD tasks and mirtask 
is used for low-level examination and manipulation of 
the data, often to accomplish tasks specific to ATA data 
reduction. 

Of particular note is that the pipeline internally man- 
ages data flow through a directed acyclic graph model 
analogous to that of the traditional Unix tool make (Feld- 
man 1979). While traditional shell scripting languages do 
not support the data structures necessary to conveniently 
implement such a design, doing so is easy in Python. The 
drPACS package (Teuben 2011) takes a similar approach 
but actually uses make for its underlying management. 
We find this approach to be unsatisfactory for two reasons. 
Firstly, make can only track dependencies between files 
on disk, and can only invoke tools via the Unix shell. 
This model maps very inefiiciently onto certain aspects of 
typical data reduction workfiows. Secondly, and more fun- 
damentally, make checks whether steps need to be rerun 
by comparing file modification times rather than actual 
contents. This suffers from a variety of minor issues and 
the major issue that if a pipeline product is regenerated, 
all products downstream of it must also be regenerated, 
even if the rebuilt product did not actually change. In a 
data reduction pipeline the costs of this inefliiciency can 
be severe. The tool described here avoids this problem 
by detecting modifications with cryptographic hashes as 
ahuded to in §2. 

7. miriad-python AS A SOFTWARE PROJECT 

In this section we describe miriad-python as a software 
project. We also hope that this section may provide 
future authors a useful reference of issues that, in our 
opinion, are important to discuss when documenting any 
other software project. 

7.1. Intellectual Property Issues 

miriad-python is an open source project, almost entirely 
licensed under the GNU General Public License, version 
3 or later. Small portions of it relating to the compilation 
process are licensed under different permissive licenses. 
As such, miriad-python is available for inspection and use 
by anyone for any purpose. The source code copyright is 
owned by the authors, miriad-python is not known to be 
subject to any software patents. 

7.2. Availability 

The source code is available for free download from the 
project website^^ either in the form of file archives or via 
the git^^ version control system. The git repository is 
also available on the website GitHub^^. Thanks to the 
distributed nature of git, clones of either repository stand 
alone and contain the complete, cryptographically- verified 
revision history of the entire project^^. Most installations 

^'^ http://purl.org/net/pkgwpub/miriad-python. This link is 
a "permanent URL" providing a durable, reroutable link suitable 
for inclusion in the academic literature. The authors encourage the 
use of this URL and not its current, possibly-ephemeral destination 
when linking to the miriad-python website and related pages. See 
http://purl.org/ for more information. 

^^ http://git-scm.com/ 

■^^ Currently https://github.com/pkgw/miriad-python/, but 
the canonical link may be found via the project website. 

^^ The SHAl checksum of the version of miriad-python described 
in this paper is aff0a3cb4d47d030426cf9e36475e4bc9ael816f. The 



of miriad-python track the git repository and so official 
releases are rare. When made, they are documented and 
linked to on the project website. The most recent release 
is version 0.6, made on 2011 April 28^^. 

Installation instructions are beyond the scope of this 
paper and are provided on the website. We note, however, 
that miriad-python must be compiled against MIRIAD 
libraries from the CARMA codebase built with the auto- 
conf-based build system. 

7.3. Development Model 

Development of miriad-python occurs in the git repos- 
itory. The authors aim to conduct development in an 
open, welcoming fashion. Contributions from the com- 
munity are encouraged via (e.g.) email or "pull requests" 
on the GitHub site. Based on the experiences of other 
community-based software projects, copyright assignment 
is not required for external contributions (see, for instance, 
the discussion in O'Mahony 2003). 

The scale of the miriad-python project is such that other 
pieces of project coordination infrastructure such as mail- 
ing lists and a bug tracker are not currently deemed nec- 
essary. These will arise as the miriad-python community 
grows, and will be linked to on the project website. 

7.4. Documentation and Examples 

A moderately-complete manual to miriad-python is avail- 
able on the website, providing both an API reference and 
somewhat higher-level guidance as to the intended usage 
of the package. The primary documentation format is 
HTML served over the web, but the documentation is gen- 
erated using Sphinx^^, the official documentation system 
of the Python language, and in principle the documen- 
tation can be generated in several other output formats, 
including printable PDF files. 

Besides the brief examples in §5, the miriad-python 
source distribution contains a small set of longer samples 
implementing features as both standalone scripts and 
importable modules. Exercised features include reading 
and writing of UV data, executing MIRIAD tasks with 
mirtask, manipulating dataset header items, reading gains 
tables, and integrating into the MIRIAD documentation 
system. 

7.5. Quality Assurance 

miriad-python is in line with the vast majority of scien- 
tific software in that, unfortunately, it has virtually no 
quality-assurance systems in place beyond the fact that it 
is exercised in day-to-day use. In particular, it has no auto- 
mated test suite. Code changes in diff format are reviewed 
before being committed to the repository, hopefully pre- 
venting unintended changes from being integrated into 
the source tree. An automatic test framework is planned, 
using short unit tests of the miriad-python APIs as well 
as more involved tests using sample MIRIAD datasets 
either generated on-the-fly (with uvgen) or optionally 
downloaded from the project website. 

secure architecture of git ensures that a vahd commit of this check- 
sum is guaranteed to correspond to the exact source code described 
in this work. 

^^ The corresponding git commit is acbfa731c5da07c0270b03e9- 
e93cllf4c9b7a4d4. 

^^ http://sphinx.pocoo.org 
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Given time constraints and the potential payoffs it is 
certainly not hard to understand why scientific software 
so often lacks systematic testing infrastructure. As a 
point of reference, testing can approach half the cost of a 
piece of software in the corporate context (Machado et al. 
2010). We emphasize, however, that software defects are 
extraordinarily common — about one for every twenty 
statement lines of code for undisciplined coding (Shull 
et al. 2002) — and that even rudimentary testing can cap- 
ture many of these defects (e.g.. Binder 1996). Shull et al. 
(2002) present an excellent summary of the prevalence 
and impact of defects in typical software systems. 

8. PERFORMANCE CONSIDERATIONS 

Python, being a dynamic interpreted language, is not 
optimal for achieving raw numerical throughput. Indeed, 
the raison d'etre of a tool such as miriad- python is that in 
many cases numerical throughput can profitably be traded 
for other desirables such as code readability, program- 
mer time savings, and ease of integration with existing 
codebases. In our experience, data-processing algorithms 
written in Python are almost always "fast enough," even 
for fairly complex algorithms. 

8.1. Serial Tests 

To give a very rough quantification of the overhead 
of using Python compared to FORTRAN-77, we wrote 
two versions of a simple task that iterates through a 
MIRIAD u-v dataset and computes the RMS value of 
the unflagged visibilities in each record. One version was 
written in Python and one in FORTRAN-77, the latter 
being compiled with GNU gfortran version 4.4.1. For a 
0.3 GB visibility dataset containing 2,812,072 u-v records 
of 16 channels each, the Python version ran about 21 
times slower than the FORTRAN-77 version (63 s vs 
3 s) on our test system, a machine running a quad-core 
AMD Opteron 1385 CPU at 2.7 GHz. (Tests were run 
repeatedly with no effort to flush filesystem caches, so 
performance was CPU-limited.) For a 1.6 GB dataset 
containing 406,260 u-v records of 1,024 channels each, the 
Python version ran about 4 times slower (26 s vs 7 s). 
These results suggest that there is indeed a nonnegligible 
overhead to the Python interpreter, although its signif- 
icance can be highly dependent on the structure of the 
input data. The performance of the Python version in the 
1,024-channel case is relatively good because relatively 
more work is done inside a few vectorized NumPy func- 
tion calls; manually iterating over each channel in Python 
increases the runtime by a factor of ~90. For cases where 
most of the work is performed inside a few vectorized 
NumPy functions, the overhead to using Python can be 
fairly small. 

In most software, the majority of the execution time is 
spent in only a small portion of the code. For those cases 
where serial throughput is a limitation, large improve- 
ments can often be achieved for low effort by identifying 
the most-executed codepaths and porting their particular 
implementations to a compiled language. This strategy 
is adopted by several of the projects discussed in §3. In 
some instances, it may be possible and profitable to gener- 
ate speed-critical code on-the-fly using a compiler library 
such as the Low-Level Virtual Machine (Lattner & Adve 
2004) , avoiding the use of customized precompiled mod- 
ules altogether. In other cases, it may make sense to 



prototype an algorithm in Python and then develop a 
production implementation in a compiled language. Of 
course, there will always exist problems for which Python 
is an inappropriate tool, but in our experience the combi- 
nation of Python and a few basic numerical libraries is 
sufficient to address a wide range of challenges. 

8.2. Parallelization 

It is clear that parallelized processing will play a fun- 
damental role in the reduction of future interferomet- 
ric datasets. There is a special opportunity for next- 
generation interferometric processing packages to estab- 
lish themselves in this space because the existing packages 
were generally built purely for serial processing. 

There are several classes of parallelization that may be 
considered. For some applications, multiple independent 
reduction processes may be executed simultaneously on a 
single machine or on a cluster with no intra-process com- 
munication or synchronization. In this "embarrassingly 
parallel" case, a language such as Python may ease the de- 
velopment of both the reduction processes and the system 
to manage the overall processing job. Indeed, there are 
several preexisting Python frameworks for constructing 
parallel applications, e.g. Parallel Python^^. 

Other parallel applications can be implemented via mul- 
tiple independent processes that communicate and syn- 
chronize. Such applications are often implemented with 
the standard Message Passing Interface (MPI; Walker 
1994), which is accessible in Python via several different 
toolkits (e.g., nnpi4py; Dalcin et al. 2008). Somewhat 
surprisingly, MPI-based parallel Python programs with 
core numerical routines written in a compiled language 
can perform as well as parallel programs written purely in 
a compiled language (Cai et al. 2005). The comparative 
ease of writing the driving logic of such a program in 
Python versus compiled languages makes this an intrigu- 
ing model for parallel algorithmic development. 

For applications that require more intensive data- 
sharing, one may wish to parallelize MIRIAD operations 
within a single process by using multithreading. Un- 
fortunately, Python's "global interpreter lock" prevents 
multiple threads from actually executing concurrently. 
(In many cases, naively multithreading a Python program 
makes it much slower!) In the particular case of miriad- 
python, extreme care would also be needed to synchronize 
access to the MIRIAD subroutine libraries, which are not 
threadsafe and maintain substantial shared state. Cotton 
(2008a) describes some of the challenges and successes 
encountered in tackling this kind of problem in Obit. The 
recommended system for threaded-style computation in 
Python is the multiprocessing module, which provides an 
infrastructure for launching concurrent Python subpro- 
cesses and communicating objects between parent and 
child. In the context of this work, multiprocessing is a 
more limited. Python-specific variant of the standard 
MPI approach. For certain I/0-limited applications, mul- 
tithreaded Python code can obtain a performance gain 
above single-threaded code, but the difficulties of seriahz- 
ing access to the MIRIAD libraries remain. 

Finally, certain classes of problems are well-suited to 
the highly-parallel capabilities of CPUs (graphics pro- 
cessing units) . Current uses of CPUs in interferometric 

^^ http : //www . parallelpython . com/ 
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appHcations include pulse dedispersion (e.g., Magro et al. 
2011; Straten & Bailes 2011), RFI mitigation (Ait-Allal 
et al. 2012), cross-correlation (Clark et al. 2011), and 
visualization (Hassan et al. 2011). Python code can take 
advantage of CPUs using frameworks such as PyCUDA 
(Klockner et al. 2012). The nature of these examples sug- 
gest that the algorithms most suited to GPU processing 
would not likely need to build on the MIRIAD/ mi riad- 
python infrastructure, except perhaps to use MIRIAD 
data formats for I/O. 

Regardless of the kind of parallelization in question, 
many interferometric algorithms are I/0-intensive. The 
options for improving the I/O bandwidth of MIRI- 
AD/miriad-python are somewhat limited due to the lack of 
built-in support for parallel I/O (e.g., Thakur et al. 1999) 
in the MIRIAD I/O subroutines. With investment in the 
appropriate hardware, speedups can be achieved using 
transparently parallel filesystems (e.g., Lustre^^ or PVFS; 
Cams et al. 2000). Significant gains in effective I/O band- 
width may alternatively be achieved with commodity 
hardware by spreading data between many independent 
processing nodes if the application permits partitioning of 
the data (e.g., by spectral channel or pointing direction). 

9. SUMMARY 

The new generation of interferometric hardware de- 
mands a new generation of interferometric software im- 
plementing a new generation of algorithms. Fortunately, 
there's a wide variety of projects aiming to provide the 
infrastructure needed to develop these tools as well as 
going ahead and implementing them. 

miriad- python is one of these projects. It is "broad" 
rather than "tall": it provides a wide range of APIs for 
accessing MIRIAD tasks and data, but does not provide 
its own algorithms built on top of those APIs. While 
certain applications will be best matched to the facilities 
provided by the higher- level packages Obit, AIPY, or 
MeqTrees, miriad-python is a good foundation on which to 
build those that are not, especially if one wishes to take 
advantage of the simple and efficient MIRIAD u-v data 
format. We actively encourage contributions to miriad- 
python development and are interested in supporting new 
miriad-python applications. 

It seems likely that SKA-scale interferometric software 
will be only tangent ially, if at all, related to any of the 
packages that exist today, and it is unclear if humans will 
"reduce" SKA data in any meaningful way given the rates 
involved. If the SKA is to meet both its goals and its 
schedule, however, the techniques needed to conquer exas- 
cale data rates must be discovered and prototyped, piece 
by piece, using today's software packages. Experience 
suggests that the innovations needed to solve these kinds 
of massive challenges do not come out of one organization; 
instead, wide-ranging and vibrant experimentation by a 
broad community of contributors will be essential. 

The authors thank W. D. Cotton for helpful discus- 
sions, miriad-python is built upon the open-source Python 
language as well as an extensive foundation of other open- 
source and Free Software packages. It could not exist 
without the work of the innumerable people whose contri- 
butions have helped build this infrastructure. Research 

^^ http://lustre.org/ 
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Foundation, the National Science Foundation, the US 
Naval Observatory, and other public and private donors. 
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Figure 1. Diagram of the relationships between different components of the miriad-python system. Applications use the three miriad-python 
modules, which in turn access the tasks and subroutines in MIRIAD, which interact with the base operating system. Both applications and 
miriad-python use the NumPy module (not depicted) for numerical array operations. 
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Figure 2. Upper panel: Spectrogram of a pulse from the Crab pulsar, summing visibilities coherently using the known location of the 
emitter. Dispersion due to the ISM causes pulse arrival times to vary with frequency (channel number). Apparent discontinuities in the 
pulse are due to flagged channels not being rendered (note the gaps in the left axis values). Lower panel: An image of the pulsar combining 
all visibilities in the observation, applying a correction for the dispersive time delay. The grayscale is linear from -42 Jy (white) to 1189 Jy 
(black). See §6.1 and Law et al. (2011). 
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Figure 3. Comparison of two methods for detecting millisecond astrophysical pulses, one using the bispectrum and one using traditional 
beamforming, both implemented in Python. Each method is applied to a visibility data stream containing five pulses from the bright pulsar 
B0329+54 and the signal-to-noise ratio (SNR) of each detection is plotted. The bispectrum-based method can be more effective (higher 
SNR detections for the same data) than coherent beamforming and is sensitive to pulses coming from any direction. See §6.1 and Law &; 
Bower (2011). 
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Figure 4. Sample screenshot of the u-v data visualizer described in §6.3. The main graphical display is a dynamic spectrum of amplitude 
on one baseline as a function of frequency (horizontal axis) and time (vertical axis). The left-hand panel provides controls for filtering and 
transforming the data. Python bindings to well-established graphical toolkits allow the data to be displayed in an attractive and responsive 
user interface. 



